#!/usr/bin/perl 

eval 'exec /usr/bin/perl  -S $0 ${1+"$@"}'
    if 0; # not running under some shell

# Copyright (C) 2009 Eric L. Wilhelm

use warnings;
use strict;

=head1 NAME

conewrap - translate geometry onto conical surface

=head1 USAGE

The cone definition is given by the apex angle (from the vertical axis
to the cone face) and the geometry is centered at a given offset from
the apex.  Output coordinates are based on the cone apex at 0,0,0.

  conewrap --angle 60d --offset 1.2  --layer 1 input.dxf > output.dxf

=head1 OPTIONS

Note that all angles are taken as radians unless given with a trailing
'd' for degrees.

=over

=item -a, --angle <angle>

Cone apex angle.  (e.g. "60d")

=item -o, --offset <distance>

Offset from apex to placement centerline.

=item -s, --arcstep <angle>

Maximum subdivision for arc approximation.  The arcs are broken into
segments (before mapping onto the cone) at or below this angle.

=item --maxphi <angle>

Maximum angle subtended by any segment around the cone.  Reduce this to
get a finer approximation of the cone surface.

=item -l, --layer <i>

Look for geometry on layer index i (default=0).  Note:  this is not the
layer name.

=item --raw

Output points in raw format.  This will be in the form:

  loop 1
    $x1,$y1,$z1
    $x2,$y2,$z2
    $x3,$y3,$z3
    ...
  loop 2
    ...

=back

=head1 See Also

L<Math::Geometry::ConeMap> has more details about the mapping method and
geometric conventions.

=cut

package bin::conewrap;

our $VERSION =
  $main::VERSION = v0.1.0;


use Getopt::Helpful;
use File::Fu;
use YAML::XS ();
use POSIX ();

use Math::Geometry::ConeMap;

use CAD::Format::DXF::Writer;
use Graph::ChainBuilder;

use constant CLASS => __PACKAGE__;

use utf8;
sub π () {atan2(1,1) * 4};

sub main {
  my (@args) = @_;

  my %o = (
    angle   => '60d',
    offset  => 1,
    onlayer => 0,
    arcstep => '15d',
    maxφ    => '5d',
  );
  my $hopt = Getopt::Helpful->new(
    usage => 'CALLER <file> [options]',
    ['a|angle=s', \$o{angle}, '<angle>',
      'apex angle (e.g. "60d")'],
    ['o|offset=f', \$o{offset}, '<distance>',
      'offset from apex to placement (centerline)'],
    ['l|layer=i',  \$o{onlayer}, '<i>',
      "use geometry on layer index i (default=$o{onlayer})"],
    ['s|arcstep=s', \$o{arcstep}, '<angle>',
      "max subdivision for arc approximation (default=$o{arcstep})"],
    ['maxphi=s', \$o{maxφ}, '<angle>',
      "maximum subtended angle (default=$o{maxφ})"],
    ['raw', \$o{raw}, '', 'output points in raw format'],
    ['version',  \$o{printv}, '', 'print version number and exit'],
    '+help',
  );
  $hopt->Get_from(\@args);
  if($o{printv}) {
    require version;
    print CLASS->VERSION, "\n";
    return;
  }

  for(qw(angle arcstep maxφ)) {
    $o{$_} *= π/180 if($o{$_} =~ s/d$//);
  }

  my $infile = shift(@args) or die "need input file";
  $infile = File::Fu->file($infile);

  my @data = grep({
    defined($_->{type}) and
    ($_->{type} eq 'arc' or $_->{type} eq 'line')
    and $_->{layer} == $o{onlayer}
  } YAML::XS::Load(do {
    my $fh = $infile->file =~ m/\.dxf/ ?
      $infile->piped_open('dxf2rzp') : $infile->open;
    local $/; readline($fh)}));

  @data or die "No entities found on layer #$o{onlayer}!\n";
  warn scalar(@data), " entities in input\n";

  my $cone = Math::Geometry::ConeMap->new(
    θ => $o{angle},
    D => $o{offset}
  );

  my $graph = Graph::ChainBuilder->new;

  for(my $i = 0; $i < @data; $i++) {
    my $ent = $data[$i];
    my $t = $ent->{type};

    flip_arc($ent) if($t eq 'arc' and $ent->{normal});

    # stringify each endpoint for optimal matching
    my @ends = CLASS->can($t . '_endpoints')->($ent);
    $_ = join(",",
      map({$_ = sprintf("%0.6f", $_); s/^-// if($_ == 0); $_} @$_[0..1])
    ) for(@ends);

    $graph->add(@ends, $ent);
  }

  if(my @ends = $graph->open_ends) {
    die "open edges found!  (dirty/duplicate entities?)",
    join("\n  ", '',
      map({join(' ', map({join("|", $_->p0, $_->p1)} @$_))} @ends)
    ), "\n";
  }

  my %functions =
    map({$_ => CLASS->can($_ . '_approximator')->($cone, %o)}
    qw(arc line));

  my %writer;
  if($o{raw}) {
    my $loop_num;
    $writer{emit} = sub {
      my @v = @_;
      $loop_num++;
      print join("\n  ", "loop $loop_num",
        map({join(",", @$_)} @v)), "\n";
    };
    $writer{finish} = sub {};
  }
  else {
    my $dw = CAD::Format::DXF::Writer->new(\*STDOUT);
    $dw->emit(0 => SECTION => 2 => ENTITIES =>);
    $writer{emit} = sub {
      my @vertices = @_;
      $dw->emit(0 => POLYLINE =>
        8 => 0, # layer
        66 => 1, # cruft
        10 => 0, 20 => 0, 30 => 0, # elevation
        70 => 9, # closed+3D flag
        map({(0 => VERTEX => 
          8 => 0, # layer
          10 => $_->[0],
          20 => $_->[1],
          30 => $_->[2],
          70 => 32, # 3D flag
        )} @vertices),
        0 => SEQEND => 8 => 0
      );
    };
    $writer{finish} = sub {
      $dw->emit(0 => ENDSEC => 0 => EOF =>);
    };
  }

  my $point_count = 0;
  foreach my $loop ($graph->loops) {
    my @vertices;
    foreach my $edge (@$loop) {
      my $ent = $edge->data;
      my $func = $functions{$ent->{type}};
      my @pts = $func->($ent);
      @pts = reverse(@pts) if($edge->reversed);
      pop(@pts); # discard redundant point
      $point_count += @pts;
      push(@vertices, @pts);
    }
    $writer{emit}->(@vertices);
  }
  $writer{finish}->();
  warn "$point_count points mapped\n";
} # end main
########################################################################

sub line_approximator {
  my $cone = shift;
  my %opts = @_;

  my $fxyz = $cone->xyz_of_xy;
  my $fZ   = $cone->z_of_xy;
  my $fφ   = $cone->φ_of_x;

  return sub {
    my $line = shift;
    my $ends = $line->{points};

    # subdivide based on dφ (from dx)
    my $dφ = $fφ->(abs($ends->[1][0] - $ends->[0][0]));
    my $n = POSIX::ceil($dφ / $opts{maxφ});

    # warn "$line->{address} dφ: ", $dφ * 180/π, "° ($n)\n"
    #   if($dφ > $opts{maxφ});

    if($n > 1) {
      @$ends = map({V(@$_)} @$ends);
      my $v = ($ends->[1] - $ends->[0]) / $n;
      foreach my $i (2..$n) {
        splice(@$ends, $i-1, 0, $ends->[0] + $v * ($i-1));
      }
    }

    my @pts = map({V($fxyz->(@$_[0,1]))} @$ends);

    # NOTE this is only check code ... can probably kill it
    # foreach my $i (1..$#pts) {
    #   my $mid = ($pts[$i] + $pts[$i-1])/2;
    #   my $diff = $fZ->(@$mid[0,1]) - $mid->[2];
    #   warn "diff: $diff\n" if($diff > 0.001);
    #   my $midmap = V($fxyz->(
    #     ($ends->[$i][0] + $ends->[$i-1][0])/2,
    #     ($ends->[$i][1] + $ends->[$i-1][1])/2,
    #   ));
    #   my $mapdiff = abs($midmap - $mid);
    #   warn "$line->{address} mapdiff: $mapdiff\n" if($mapdiff > 0.002);
    # }

    return @pts;
  };
}
########################################################################

sub arc_approximator {
  my $cone = shift;
  my %opt = @_;

  my $fxyz = $cone->xyz_of_xy;
  my $fZ   = $cone->z_of_xy;

  return sub {
    my $arc = shift;

    my $r = $arc->{radius};
    my $c = V(@{$arc->{points}[0]});
    my ($a0, $a1) = @{$arc->{angles}};

    my $delta = $a1 - $a0;
    $delta = π * 2 + $delta if($delta < 0);

    my @pts;
    # start by breaking at 15° or less?
    my $n = POSIX::ceil($delta / ($opt{arcstep}));
    my $step = $delta / $n;
    foreach my $i (0..$n) {
      my $A = $a0 + $step * $i;
      my $pt = $c + V($r * cos($A), $r * sin($A));
      # TODO check the 2D dx as dφ vs maxφ here and sub-segment?
      #  (or do it in another pass later)
      $pt = V($fxyz->(@$pt[0,1]));
      push(@pts, $pt);
    }

    # TODO we might need to actually apply the maxphi limit here
    # this is just sanity checking
    foreach my $i (1..$#pts) {
      my $mid = ($pts[$i] + $pts[$i-1])/2;
      my $diff = $fZ->(@$mid[0,1]) - $mid->[2];
      warn "depth: $diff\n" if($diff > 0.001);
      my $dφ = abs(
        atan2($pts[$i][1], $pts[$i][0]) -
        atan2($pts[$i-1][1], $pts[$i-1][0])
      );
      $dφ -= π if($dφ > π);
      warn "arc $arc->{address} ($i) ", $dφ*180/π, "°\n" if($dφ > $opt{maxφ});
    }

    return @pts;
  };
}
########################################################################

use Math::Vec qw(V);

# only dealing with flipping angles across the y-axis here
sub flip_arc {
  my $ent = shift;

  my $norm = delete($ent->{normal});
  croak("unexpected @$norm") unless(join(',', @$norm) eq '0.0,0.0,-1.0');
  foreach my $ang (@{$ent->{angles}}) {
    $ang = $ang > π ? 3*π-$ang : π-$ang;
  }
  @{$ent->{angles}} = reverse @{$ent->{angles}};
  return($ent);
}
########################################################################

sub arc_endpoints {
  my ($arc) = @_;

  my $r = $arc->{radius};
  my $c = V(@{$arc->{points}[0]});
  #warn "c: $c";
  my ($a0, $a1) = @{$arc->{angles}};
  return map({$c + V($r * cos($_), $r * sin($_))} $a0, $a1);
}
########################################################################

sub line_endpoints {
  my ($line) = @_;

  return map({V(@$_)} @{$line->{points}});
}
########################################################################

package main;

if($0 eq __FILE__) {
  bin::conewrap::main(@ARGV);
}

# vi:ts=2:sw=2:et:sta
my $package = 'bin::conewrap';
